#include "head.h"
#include "vect_3d.cpp"
#include "class_ma.cpp"
COMPLEX _i=COMPLEX(0,1);
double pi=3.1415926535897932384626433832795;
                                                                  
//The following variables are used in mgs().
//mgs==modified Gram Schmidt method for orthogonal factorization.
const int largest_array_size=10000;
matrix_ sk=matrix_(largest_array_size,1),uk=matrix_(largest_array_size,1);
double norm_a[largest_array_size];
Complex vv[largest_array_size],uu[largest_array_size],ss[largest_array_size],rr[largest_array_size],qq[largest_array_size];

//This is the modified MGS2: All the arrays are substituted by the pointers.
int MGS2(matrix_ &AA,matrix_ &QQ,matrix_ &RR,double er,bool &switcher)
{
	int m=AA.rows,n=AA.columns;
   int k,l,j,k1;
   double max_norm,f_norm0,f_norm1,err=er*er;
   int max_norm_k,kn;
   int r;
   int min_2=(m>n)?n/2:m/2;
   Complex ss1;

   double ETA=0.1;  //0.5; //0.707*0.707;
   double ratio,ratio0;
   COMPLEX temp;
   COMPLEX *A=AA.array,*Q=QQ.array,*R=RR.array,*t,*t1,*t2,*t3;
   double *tt;
   Complex *t4;

   copy(QQ,0);copy(RR,0);

   //Notice that any norm here is the square of the norm defined in math book.
   //We use this definition just for the sake of convenience.

   max_norm=-1;max_norm_k=0;f_norm0=0;
   for(k=0;k<n;k++)                                    //
   {                                                   //Compute the coulmn norms
      temp=COMPLEX(0);                                 //
      t=A+k;			                                 //of A.
   	for(j=0;j<m;j++)
      {  					                               //
	      temp+=conj(*t)**t;
         t+=n;
      }
      tt=norm_a+k;										         //and the Frobenious Norm
      *tt=abs(temp);                           			//of A.
      f_norm0+=*tt;

      if(max_norm<*tt)                                 //Seek the column having
      {                                                //the maximun norm in A.
      	max_norm=*tt;                                 //Store its column number
      	max_norm_k=k;                                 //into max_norm_k.
      }
   }
   if(max_norm<1e-20)
   {
   	switcher==true;
      return 1;
   }


   t=A;
   t1=A+max_norm_k;
   for(j=0;j<m;j++)                                    // Interchange column
   {                                                   // max_norm_k with coulmn
      temp=*t;		                                     // 0 of A.
      *t=*t1; 						                       	 //
      *t1=temp;
      t+=n;t1+=n;
   }                          								 //

   *(norm_a+max_norm_k)=*norm_a;                    	 // Simultaneously,
   *norm_a=max_norm;		                               // Interchange their norms.


   r=0;kn=0;
   for(k=0;k<n;k++)
   {
      t=R+kn+k;t1=A+k;
      for(l=0;l<m;l++)                                //  Calculate the
      {
   		*t+=*t1*conj(*t1);   								//  diagonal entries
         t1+=n;
      }
      *t=sqrt(*t);				                        //  of matrix R.


      t1=A+k;t2=Q+k;
      for(l=0;l<m;l++)                                //  Comput the k th
      {
      	*t2=*t1/(*t);				                     //  column of matrix Q.
         t1+=n;t2+=n;
      }
      r+=1;

 


      k1=k+1;max_norm=-1;max_norm_k=0; f_norm1=0;

      for(j=k1;j<n;j++)
      {
         t1=R+kn+j;t2=Q+k;t3=A+j;
         for(l=0;l<m;l++)                             //  Compute the entry
         {
      		*t1+=conj(*t2)**t3;   							//  R[k][j] of matrix R.
            t2+=n;t3+=n;
         }

         t2=A+j;t3=Q+k;
         for(l=0;l<m;l++)                             //  Overwrite the j th
         {
         	*t2-=*t3**t1;       								//  column of A by substrating
            t2+=n;t3+=n;
                                                      //  its partial along the direction
                                                      //  of the k th column of Q.
         }

         t2=A+j;t3=Q+k;
         ratio=0;
         for(l=0;l<m;l++)
         {
         	ratio=ratio+abs(conj(*t2)**t2);
            *(qq+l)=*t3;
            t2+=n;t3+=n;
         }

         if(ratio<ETA**(norm_a+j))
         {
            t2=A+j;
         	for(l=0;l<m;l++)
            {
            	*(vv+l)=*t2;
               t2+=n;
            }

            do
            {
               for(l=0;l<m;l++)*(uu+l)=Complex(0);

            	ss1=Complex(0);
               for(l=0;l<m;l++)
               ss1=ss1+conj(*(qq+l))**(vv+l);        //s=qH*aj;

               for(l=0;l<m;l++)
               *(uu+l)+=*(qq+l)*ss1;

               ratio0=ratio;
               ratio=0;

               for(l=0;l<m;l++)
               {
                  t4=vv+l;
               	*t4-=*(uu+l);                        //v=v-u;
                  ratio+=abs(conj(*t4)**t4);
               }
            }
            while(ratio<ETA*ratio0);

            t2=A+j;
            for(l=0;l<m;l++)
            {
            	*t2=*(vv+l);
               t2+=n;
            }
         }

         tt=norm_a+j;
         *tt=ratio;

         f_norm1+=*tt;
         if(max_norm<*tt)
         {
         	max_norm=*tt;
            max_norm_k=j;
         }
      }

      if(f_norm1>err*f_norm0)//&&r<min_2)
      {
         t3=A+k1;t2=A+max_norm_k;
      	for(l=0;l<m;l++)
      	{
      		temp=*t3;
         	*t3=*t2;
         	*t2=temp;
            t2+=n;t3+=n;
      	}
         *(norm_a+max_norm_k)=*norm_a;
         *(norm_a+k1)=max_norm;

      }
      else break;
      kn+=n;
   }

   if(r>=min_2){switcher=false;}
   else
   {
   	switcher=true;
   }
   return r;
}


//////////////////////////////////////////////////////////////////////////////////////
//Function 11. lagrange interpolation base function.
//j: the base index.
//n: the number of zeros.
//tk: the array stroing the zeros.
//t: the point that should be interpolated.
float lagb(int j,int n,float *tk,float t) //lagrange base function.
{
   int k;
   float f;
   f=1;
	for(k=0;k<n;k++)
   {
   	if(j!=k)f*=(t-tk[k])/(tk[j]-tk[k]);
   }
   return f;
}

////////////////////////////////////////////////////////////////////////////////////
//Function 12. Lagrange base function of three dimensional problem.
//n: number of zeros in one dimension.
//i,j,k: the indices of the weight.
//xi,yi,zi the coordinates of the zeros.
//x,y,z: the point should be interpolated.
float lagbase3d(int i,int j,int k,float *xi,float *yi,float *zi,float x,float y,float z,int n)
{
	float f;
   f=lagb(i,n,xi,x)*lagb(j,n,yi,y)*lagb(k,n,zi,z);
   return f;
}

float sinc(float tj,float t,float pi_dt)
{
   float f;
	float x=(t-tj)*pi_dt;
   if(fabs(x)<1e-20)f=1;
   else f=sin(x)/x;
   return f;
}
float sincd(float tj,float t,float pi_dt,int N)
{
   float f;
	float x=(t-tj)*pi_dt;
   if(fabs(x)<1e-20)f=1;
   else f=sin(x)/N/sin(x/N);
   return f;
}
float sinc3d(float xi,float yi,float zi,float x,float y,float z,float pi_dx,float pi_dy,float pi_dz)
{
	float f;
   f=sinc(xi,x,pi_dx)*sinc(yi,y,pi_dy)*sinc(zi,z,pi_dz);
   return f;
}

float sincd3d(float xi,float yi,float zi,float x,float y,float z,
float pi_dx,float pi_dy,float pi_dz,int Nx,int Ny,int Nz)
{
	float f;
   f=sincd(xi,x,pi_dx,Nx)*sincd(yi,y,pi_dy,Ny)*sincd(zi,z,pi_dz,Nz);
   return f;
}

void Green(matrix_ &G,vector &dir,float dl, int K1,float k)
{
   int n,o,p,q,r,s,u;
   vector rf,rs;
   double R;

   u=0;
   for(n=0;n<K1;n++)
   for(o=0;o<K1;o++)
   for(p=0;p<K1;p++)
   {
   	rf=dir+vector(n,o,p)*dl;
   	for(q=0;q<K1;q++)
   	for(r=0;r<K1;r++)
   	for(s=0;s<K1;s++)
   	{
   		rs=vector(q,r,s)*dl;
         R=(rf-rs).abs_vector();
   		G[u]=exp(_i*k*R)/R; //The type of green function can be changed.
   		u++;
      }
   }
}

void Green(matrix_ &G,float dl, int K1,float k,vector rs)
{
   int n,o,p,q,r,s,u;
   vector rf;
   double R;

   u=0;
   for(n=0;n<K1;n++)
   for(o=0;o<K1;o++)
   for(p=0;p<K1;p++)
   {
   	rf=vector(n,o,p)*dl;
      R=(rf-rs).abs_vector();
      G[u]=exp(_i*k*R)/R; //The type of green function can be changed.
      u++;
   }
}

void WWW(matrix_ &W,vector *r,int N,float *xi,float *yi,float *zi,int K1)
{
	int n,i,j,k,u;
   u=0;
   for(n=0;n<N;n++)
   {
   	for(i=0;i<K1;i++)
   	for(j=0;j<K1;j++)
   	for(k=0;k<K1;k++)
   	{
      	W[u]=COMPLEX(lagbase3d(i,j,k,xi,yi,zi,r[n].vect_x,r[n].vect_y,r[n].vect_z,K1));          
         u++;
      }
   }
}
void WWW_sinc(matrix_ &W,vector *r,int N,float *xi,float *yi,float *zi,int K1)
{
	int n,i,j,k,u;
   float xm,ym,zm,pi_dx,pi_dy,pi_dz;
   pi_dx=pi/(xi[1]-xi[0]);
   pi_dy=pi/(yi[1]-yi[0]);
   pi_dz=pi/(zi[1]-zi[0]);
   u=0;
   for(n=0;n<N;n++)
   {
   	for(i=0;i<K1;i++)
   	for(j=0;j<K1;j++)
   	for(k=0;k<K1;k++)
   	{
         xm=xi[i];ym=yi[j];zm=zi[k];
      	W[u]=COMPLEX(sinc3d(xm,ym,zm,r[n].vect_x,r[n].vect_y,r[n].vect_z,pi_dx,pi_dy,pi_dz));
         u++;
      }
   }
}

void WWW_sincd(matrix_ &W,vector *r,int N,float *xi,float *yi,float *zi,int K1)
{
	int n,i,j,k,u;
   float xm,ym,zm,pi_dx,pi_dy,pi_dz;
   pi_dx=pi/(xi[1]-xi[0]);
   pi_dy=pi/(yi[1]-yi[0]);
   pi_dz=pi/(zi[1]-zi[0]);
   u=0;
   for(n=0;n<N;n++)
   {
   	for(i=0;i<K1;i++)
   	for(j=0;j<K1;j++)
   	for(k=0;k<K1;k++)
   	{
         xm=xi[i];ym=yi[j];zm=zi[k];
      	W[u]=COMPLEX(sincd3d(xm,ym,zm,r[n].vect_x,r[n].vect_y,r[n].vect_z,pi_dx,pi_dy,pi_dz,K1,K1,K1));
         u++;
      }
   }
}

void test1() //the source points are not fixed
{
   float length=1,lambda=0.5,kin=2*pi/lambda;
	vector dir=vector(2*length,0,0);
   int K1=3,N1=6,K3D=K1*K1*K1,N3D=N1*N1*N1;
   float dl=length/(K1-1),dln=length/(N1-1);

   float *xi,*yi,*zi;
   vector *rr;
   rr=new vector[N3D];
   xi=new float[K1];
   yi=new float[K1];
   zi=new float[K1];

   matrix_ Wf,Ws,GG;

   Wf=matrix_(N3D,K3D);
   Ws=matrix_(N3D,K3D);
   GG=matrix_(K3D,K3D);

   int i,j,k,o,p,q,u,v,r,s;

   //generate the source points.
   u=0;
   for(i=0;i<N1;i++)
   for(j=0;j<N1;j++)
   for(k=0;k<N1;k++)
   {
   	rr[u]=vector(i*dln,j*dln,k*dln);
      u++;
   }

   //generate the source interpolation points.
   for(i=0;i<K1;i++)
   {
   	xi[i]=yi[i]=zi[i]=i*dl;
   }

   //generate the interpolation matrix for the source region.
   WWW(Ws,rr,N3D,xi,yi,zi,K1);
   //WWW_sinc(Ws,rr,N3D,xi,yi,zi,K1);
   //WWW_sincd(Ws,rr,N3D,xi,yi,zi,K1);
   //generate the field points.
   u=0;
   for(i=0;i<N1;i++)
   for(j=0;j<N1;j++)
   for(k=0;k<N1;k++)
   {
   	rr[u]=dir+vector(i*dln,j*dln,k*dln);
      u++;
   }

   //generate the field interpolation points.
   for(i=0;i<K1;i++)
   {
   	xi[i]=dir.vect_x+i*dl;
      yi[i]=dir.vect_y+i*dl;
      zi[i]=dir.vect_z+i*dl;
   }

   //generate the interpolation matrix for the feild region.
   WWW(Wf,rr,N3D,xi,yi,zi,K1);
   //WWW_sincd(Wf,rr,N3D,xi,yi,zi,K1);

   //generate the green's function.
   Green(GG,dir,dl,K1,kin);


   /*
   matrix_ GG1,GG2,dG;
   GG1=matrix_(K3D,K3D);
   GG2=matrix_(K3D,K3D);
   copy(GG1,GG);
   H(GG,GG2);

   matrix_ QQ,RR,QQ1,RR1,QQ2,RR2;
   QQ=matrix_(K3D,K3D);RR=matrix_(K3D,K3D);
   QQ1=matrix_(K3D,K3D);RR1=matrix_(K3D,K3D);
   QQ2=matrix_(K3D,K3D);RR2=matrix_(K3D,K3D);

   int rank1,rank2;
   bool switcher=false;

   rank1=MGS2(GG1,QQ1,RR1,1e-3,switcher);
   switcher=false;
   rank2=MGS2(GG2,QQ2,RR2,1e-3,switcher);
   compression(QQ1,QQ1.rows,rank1);
   compression(QQ2,QQ2.rows,rank2);

   change_rows_cols(RR1,rank1,K3D);
   change_rows_cols(RR2,rank2,K3D);

   H(GG,GG2);
   Hmulti(RR1,QQ1,GG);
   Hmulti(RR2,QQ2,GG2);


   matrix_ SS,SSH;

   SSH=matrix_(rank2,rank1);
   multi(SSH,RR2,QQ1);
   change_rows_cols(GG1,QQ2.rows,rank1);
   multi(GG1,QQ2,SSH);
   SS=matrix_(GG1.columns,GG1.rows);
   H(GG1,SS);
   multi(GG2,QQ1,SS);

   double er;
   dG=matrix_(K3D,K3D);
   sub(dG,GG,GG2);
   er=abs_matrix(dG)/abs_matrix(GG);
   cout<<"\nerror="<<er;

   change_rows_cols(QQ1,rank2,rank1);
   change_rows_cols(RR1,rank1,rank1);
   rank2=MGS2(SSH,QQ1,RR1,1e-3,switcher);
   cout<<"\nrank="<<rank2;
   */





   //test if the interpolation is accurate.
   double errormax=-10000,errormin=1e27, erroravg=0,error,R;
   COMPLEX G0,G;
   vector rs,rf;

   u=0;
   for(i=0;i<N1;i++)
   for(j=0;j<N1;j++)
   for(k=0;k<N1;k++)
   {
      v=0;
      rf=dir+vector(i*dln,j*dln,k*dln);
      //cout<<"\ni,j,k="<<i<<' '<<j<<' '<<k<<' ';
   	//for(o=0;o<N1;o++)
      //for(p=0;p<N1;p++)
     // for(q=0;q<N1;q++)
      {
         o=N1-1;p=N1/2;q=N1/2;
         v=o*N1*N1+p*N1+q;
      	rs=vector(o*dln,p*dln,q*dln);
         //rs=vector(length,0.5*length,0.5*length);

         R=(rf-rs).abs_vector();
         G0=exp(_i*kin*R)/R;
         G=COMPLEX(0);
         for(r=0;r<K3D;r++)
         for(s=0;s<K3D;s++)
         {
         	G=G+Wf[u*K3D+r]*GG[r*K3D+s]*Ws[v*K3D+s];
         }

         error=abs((G-G0)/G0);
         //cout<<"\nerror="<<error;
         if(error>errormax)errormax=error;
         if(errormax>1)
         {
         cout<<"stop";
         }
         if(error<errormin)errormin=error;
         erroravg+=error;
         v++;
      }
      u++;
   }
   erroravg/=N3D;//*N3D;
   cout<<"\nerrormax="<<errormax;
   cout<<"\nerrormin="<<errormin;
   cout<<"\nerroravg="<<erroravg;
   char key;
   cin>>key;
}

//the interaction of A point to a cube.
//test Lagrange interpolation with tartan grid
void test2()
{
   float length=2,lambda=1,kin=2*pi/lambda;
	vector dir=vector(2*length,0,0);
   int K1=3,N1=6,K3D=K1*K1*K1,N3D=N1*N1*N1;
   float dl=length/(K1-1),dln=length/(N1-1);
   ofstream fout;
   fout.open("error_lag.txt");

   for(K1=3;K1<=10;K1++)
   {
   N1=20;
   K3D=K1*K1*K1;
   N3D=N1*N1*N1;
   dl=length/(K1-1);dln=length/(N1-1);
   float *xi,*yi,*zi;
   vector *rr;
   rr=new vector[N3D];
   xi=new float[K1];
   yi=new float[K1];
   zi=new float[K1];

   matrix_ Wf,GG;

   Wf=matrix_(N3D,K3D);
   GG=matrix_(K3D,1);

   int i,j,k,o,p,q,u,v,r,s;


   //generate the field points.
   u=0;
   for(i=0;i<N1;i++)
   for(j=0;j<N1;j++)
   for(k=0;k<N1;k++)
   {
   	rr[u]=vector(i*dln,j*dln,k*dln);
      u++;
   }

   //generate the field interpolation points.
   for(i=0;i<K1;i++)
   {
   	xi[i]=i*dl;
      yi[i]=i*dl;
      zi[i]=i*dl;
   }

   //generate the interpolation matrix for the field region.
   WWW(Wf,rr,N3D,xi,yi,zi,K1);
   //WWW_sincd(Wf,rr,N3D,xi,yi,zi,K1);

   vector rs;
   rs=dir+vector(0,0.5*length,0.5*length);
   //generate the green's function.
   Green(GG,dl,K1,kin,rs);


   //test if the interpolation is accurate.
   double errormax=-10000,errormin=1e27, erroravg=0,error,R;
   COMPLEX G0,G;
   vector rf;

   u=0;
   for(i=0;i<N1;i++)
   for(j=0;j<N1;j++)
   for(k=0;k<N1;k++)
   {
      rf=vector(i*dln,j*dln,k*dln);
      {

         R=(rf-rs).abs_vector();
         G0=exp(_i*kin*R)/R;
         G=COMPLEX(0);
         for(r=0;r<K3D;r++)
         {
         	G=G+Wf[u*K3D+r]*GG[r];
         }

         error=abs((G-G0)/G0);
         //cout<<"\nerror="<<error;
         if(error>errormax)errormax=error;
         if(error<errormin)errormin=error;
         erroravg+=error;
         v++;
      }
      u++;
   }
   erroravg/=N3D;
   cout<<"\nerroravg="<<erroravg;
   fout<<K3D<<' '<<erroravg<<'\n';
   }
}

COMPLEX expansion_of_G0(float k,float R,int N)
{
	COMPLEX G0,i=COMPLEX(0,1),term,ikR,ikR_n;
   int n;
   float n_factorial;
   G0=term=COMPLEX(1);
   ikR=i*k*R;
   for(n=1;n<=N;n++)
   {
      ikR_n=ikR/float(n);
      term=term*ikR_n;
      G0=G0+term;
   }
   G0=G0/R;
   return G0;

}

void test4()
{
	float R,k,lambda,pi=3.1415926535897932384626433832795;
   int n,N=50;
   lambda=1;k=2*pi/lambda;R=3;
   COMPLEX G0,G,i=COMPLEX(0,1);
   G=expansion_of_G0(k,R,N);
   G0=exp(i*k*R)/R;
   cout<<"G="<<G<<" G0="<<G0;
   char key;
   cin>>key;
}
void main()
{
	//test1(); //source region is a cube
   test2();   //source region is a point and use lagrange interpolation

}




